#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Mesh generation script for the Dragonfly Aerodynamics Simulation.

This script uses GMSH to create an idealized dragonfly wing mesh based on 
the specifications in the report:
- Zero-thickness flat quarter-ellipse planform
- Span: 7 cm (0.07 m)
- Root chord: 4 cm (0.04 m)
- Hybrid mesh with structured zone around wing and unstructured outer zone
"""

import gmsh
import numpy as np
import os
import sys

# Add the project root to Python path
PROJECT_ROOT = os.path.dirname(os.path.abspath(__file__))
sys.path.append(PROJECT_ROOT)

import config

def create_dragonfly_mesh(output_file="mesh/dragonfly_hybrid_mesh.msh"):
    """
    Creates a mesh for the dragonfly wing simulation using GMSH.
    
    Args:
        output_file (str): Path to save the mesh file
    """
    # Initialize GMSH
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 1)
    gmsh.model.add("dragonfly_wing")
    
    # Get geometry parameters from config
    span = config.WING_SPAN  # 0.07 m
    root_chord = config.WING_ROOT_CHORD  # 0.04 m
    
    # Domain dimensions (much larger than wing)
    domain_length = 20 * root_chord  # 0.8 m
    domain_width = 10 * span  # 0.7 m
    domain_height = 10 * span  # 0.7 m
    
    # --- Create the wing geometry ---
    # Quarter-ellipse for the wing planform (zero thickness)
    # We'll create it in the x-y plane centered at origin
    
    # Points for the quarter-ellipse
    # Using parametric definition: x = a*cos(t), y = b*sin(t)
    # where a = span/2, b = root_chord/2, t from 0 to pi/2
    
    n_wing_points = 20  # Number of points on wing outline
    wing_points = []
    
    # Create wing outline points
    for i in range(n_wing_points):
        t = i * (np.pi/2) / (n_wing_points - 1)
        x = (span/2) * np.cos(t)
        y = (root_chord/2) * np.sin(t) 
        z = 0.0
        point_tag = gmsh.model.geo.addPoint(x, y, z, 0.002)  # Small mesh size on wing
        wing_points.append(point_tag)
    
    # Add root point (0, 0, 0)
    root_point = gmsh.model.geo.addPoint(0, 0, 0, 0.002)
    
    # Create wing outline curves
    wing_curves = []
    
    # Curve along the ellipse
    for i in range(len(wing_points) - 1):
        curve = gmsh.model.geo.addLine(wing_points[i], wing_points[i+1])
        wing_curves.append(curve)
    
    # Closing curves
    curve_tip_to_root = gmsh.model.geo.addLine(wing_points[-1], root_point)
    curve_root_to_leading = gmsh.model.geo.addLine(root_point, wing_points[0])
    wing_curves.extend([curve_tip_to_root, curve_root_to_leading])
    
    # Create wing surface (zero thickness)
    wing_loop = gmsh.model.geo.addCurveLoop(wing_curves)
    wing_surface = gmsh.model.geo.addPlaneSurface([wing_loop])
    
    # --- Create the computational domain ---
    # Box domain centered at origin
    box_points = [
        gmsh.model.geo.addPoint(-domain_length/2, -domain_width/2, -domain_height/2, 0.1),
        gmsh.model.geo.addPoint( domain_length/2, -domain_width/2, -domain_height/2, 0.1),
        gmsh.model.geo.addPoint( domain_length/2,  domain_width/2, -domain_height/2, 0.1),
        gmsh.model.geo.addPoint(-domain_length/2,  domain_width/2, -domain_height/2, 0.1),
        gmsh.model.geo.addPoint(-domain_length/2, -domain_width/2,  domain_height/2, 0.1),
        gmsh.model.geo.addPoint( domain_length/2, -domain_width/2,  domain_height/2, 0.1),
        gmsh.model.geo.addPoint( domain_length/2,  domain_width/2,  domain_height/2, 0.1),
        gmsh.model.geo.addPoint(-domain_length/2,  domain_width/2,  domain_height/2, 0.1)
    ]
    
    # Create box edges
    box_lines = [
        # Bottom face
        gmsh.model.geo.addLine(box_points[0], box_points[1]),
        gmsh.model.geo.addLine(box_points[1], box_points[2]),
        gmsh.model.geo.addLine(box_points[2], box_points[3]),
        gmsh.model.geo.addLine(box_points[3], box_points[0]),
        # Top face
        gmsh.model.geo.addLine(box_points[4], box_points[5]),
        gmsh.model.geo.addLine(box_points[5], box_points[6]),
        gmsh.model.geo.addLine(box_points[6], box_points[7]),
        gmsh.model.geo.addLine(box_points[7], box_points[4]),
        # Vertical edges
        gmsh.model.geo.addLine(box_points[0], box_points[4]),
        gmsh.model.geo.addLine(box_points[1], box_points[5]),
        gmsh.model.geo.addLine(box_points[2], box_points[6]),
        gmsh.model.geo.addLine(box_points[3], box_points[7])
    ]
    
    # Create box faces
    box_faces = []
    # Bottom face
    bottom_loop = gmsh.model.geo.addCurveLoop([box_lines[0], box_lines[1], box_lines[2], box_lines[3]])
    box_faces.append(gmsh.model.geo.addPlaneSurface([bottom_loop]))
    # Top face
    top_loop = gmsh.model.geo.addCurveLoop([box_lines[4], box_lines[5], box_lines[6], box_lines[7]])
    box_faces.append(gmsh.model.geo.addPlaneSurface([top_loop]))
    # Side faces
    for i in range(4):
        j = (i + 1) % 4
        side_loop = gmsh.model.geo.addCurveLoop([box_lines[i], box_lines[8+j], -box_lines[4+i], -box_lines[8+i]])
        box_faces.append(gmsh.model.geo.addPlaneSurface([side_loop]))
    
    # Create box volume
    surface_loop = gmsh.model.geo.addSurfaceLoop(box_faces)
    domain_volume = gmsh.model.geo.addVolume([surface_loop])
    
    # --- Define physical groups (tags) ---
    # These are crucial for the simulation to identify different parts of the mesh
    
    # Wing surface
    gmsh.model.addPhysicalGroup(2, [wing_surface], name="wing_surface")
    
    # Domain boundaries
    inlet_face = box_faces[4]  # -x face (assuming flow in +x direction)
    outlet_face = box_faces[5]  # +x face
    gmsh.model.addPhysicalGroup(2, [inlet_face], name="inlet")
    gmsh.model.addPhysicalGroup(2, [outlet_face], name="outlet")
    gmsh.model.addPhysicalGroup(2, box_faces[0:4] + box_faces[6:], name="farfield")
    
    # Volume
    gmsh.model.addPhysicalGroup(3, [domain_volume], name="fluid_domain")
    
    # --- Create structured zone around wing ---
    # For simplicity, we'll create a cylindrical region around the wing
    # that will have finer mesh
    inner_radius = 3 * root_chord  # 0.12 m
    inner_cylinder_points = []
    n_cylinder = 16
    
    # Create cylinder points at z=0
    for i in range(n_cylinder):
        theta = 2 * np.pi * i / n_cylinder
        x = inner_radius * np.cos(theta)
        y = inner_radius * np.sin(theta)
        z = 0.0
        point = gmsh.model.geo.addPoint(x, y, z, 0.01)  # Finer mesh in inner zone
        inner_cylinder_points.append(point)
    
    # Create cylinder curves
    inner_cylinder_curves = []
    for i in range(n_cylinder):
        j = (i + 1) % n_cylinder
        curve = gmsh.model.geo.addLine(inner_cylinder_points[i], inner_cylinder_points[j])
        inner_cylinder_curves.append(curve)
    
    # --- Synchronize before meshing ---
    gmsh.model.geo.synchronize()
    
    # --- Set mesh refinement fields ---
    # Distance field from wing
    gmsh.model.mesh.field.add("Distance", 1)
    gmsh.model.mesh.field.setNumbers(1, "CurvesList", wing_curves)
    gmsh.model.mesh.field.setNumber(1, "Sampling", 100)
    
    # Threshold field for mesh size based on distance
    gmsh.model.mesh.field.add("Threshold", 2)
    gmsh.model.mesh.field.setNumber(2, "InField", 1)
    gmsh.model.mesh.field.setNumber(2, "SizeMin", 0.002)  # Min size near wing
    gmsh.model.mesh.field.setNumber(2, "SizeMax", 0.05)   # Max size far from wing
    gmsh.model.mesh.field.setNumber(2, "DistMin", 0.01)   # Distance where min size starts
    gmsh.model.mesh.field.setNumber(2, "DistMax", 0.2)    # Distance where max size is reached
    
    gmsh.model.mesh.field.setAsBackgroundMesh(2)
    
    # --- Generate the mesh ---
    print("Generating 3D mesh...")
    gmsh.model.mesh.generate(3)
    
    # --- Save the mesh ---
    output_dir = os.path.dirname(output_file)
    if output_dir and not os.path.exists(output_dir):
        os.makedirs(output_dir)
    
    gmsh.write(output_file)
    print(f"Mesh saved to: {output_file}")
    
    # Get some statistics
    nodes = gmsh.model.mesh.getNodes()
    print(f"Total nodes: {len(nodes[0])}")
    
    # Finalize GMSH
    gmsh.finalize()


if __name__ == "__main__":
    # Check if gmsh is installed
    try:
        import gmsh
    except ImportError:
        print("ERROR: GMSH Python API is not installed.")
        print("Please install it using: pip install gmsh")
        sys.exit(1)
    
    # Generate the mesh
    create_dragonfly_mesh() 